Non different iable large-deviation functionals in boundary-driven diffusive systems 
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We study the probability of arbitrary density profiles in conserving diffusive fields which are driven 
by the boundaries. We demonstrate the existence of singularities in the large-deviation functional, 
the direct analog of the free-energy in non-equilibrium systems. These singularities are unique to 
non-equilibrium systems and are a direct consequence of the breaking of time-reversal symmetry. 
This is demonstrated in an exactly-solvable model and also in numerical simulations on a boundary- 
driven Ising model. We argue that this singular behavior is expected to occur in models where the 
compressibility has a deep enough minimum. The mechanism is explained using a simple model. 
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Consider a field p(x) with diffusive dynamics which 
are conserving in the bulk. Here, p(x) could describe the 
density of a gas in a capillary connecting two reservoirs. 
When the reservoirs are of equal density p, the system is 
in equilibrium. Then the steady-state probability P[pf] 
of an arbitrary density profile Pf(x) is given by -P[p/] ~ 
e - p {ps\/ k BT ^ w h ere p j s tfi e free-energy. Generically, for 
a system with short-range interactions, F is a local func- 
tional of Pf(x) and, in the disordered phase, it is a smooth 
functional. For instance, when the particles in the cap- 
illary interact only through hard-core exclusion, F[pf\ = 

NJo dx {p f (x) log + Pf (x)) log i^M}, where 
N is the length of the capillary and the density is nor- 
malized such that p = 1 corresponds to a filled capillary. 

A fundamental goal of non-equilibrium statistical me- 
chanics is to evaluate and understand the general struc- 
ture of P[pf] when the densities of the two reservoirs are 
different, such that a current flows through the system. 
For diffusive systems there has been great progress in re- 
cent years, and by now several important properties of 
P\fif] have been established. For example, it is known 
that P[pf] ~ e~ Nd ^ p f\ where N d is the volume of the 
system and <f>\pf] is a non-equilibrium analogue of the 
free energy, called the large deviation functional (LDF). 
It attains a minimum at the most probable density pro- 
file and, when the system is driven out of equilibrium, it 
becomes a non-local functional of pf (x). This has impor- 
tant consequences, manifested for example through long 
range correlations which arepresent even when the inter- 
actions are strictly local [H,0- I n contrast to equilibrium 
the LDF depends on the dynamics, and not only on the 
Boltzmann weights. 

Recently, framework within which <f>[pf] can be calcu- 
lated has been laid out, building on standard tools from 
the theory of large deviations [3|-|8[ • The key observation 
is that in these systems the probability of an atypical 
event is dominated by a single history leading up to it, 
and starting from the most probable profile; other histo- 
ries are exponentially less likely in the system size. While 
calculating the LDF within the framework is in general 
extremely difficult, it has allowed to establish exact solu- 
tions in a number of simple models [3, |j| (in some cases 
building on solutions obtained by other methods 10]). 



In addition, the framework has led to efficient numeri- 
cal algorithms for evaluating P[pf] for arbitrary diffusive 
models [ll|, as well as LDFs of global quantities such as 
the current (H-fijl. 



Despite the successes, a general understanding of the 
properties of P[pf] out of equilibrium is lacking. For 
example, for boundary-driven currents induced by reser- 
voirs held at different densities, 4>[pf] is a smooth func- 
tional of pj for particles diffusing with hard-core exclu- 
sion, just as in equilibrium 110]. In stark contrast, in 
the presence of a strong bulk drive (for example, when a 
constant force acts on the particles in the capillary) <fi[pf] 
becomes non-analytic [l6| . This is due to the existence of 
multiple histories leading to the same pf with compara- 
ble weight, and immediately implies singular behavior in 
the LDF of other more global quantities. These singular- 
ities are more subtle than those observed for the current 
or particle number [l7l - |20l ]. which can arise even if the 
LDF on the phase space is smooth, hence even in equi- 
librium models (in contrast to the present phenomenon, 
see below). It is far from clear which models exhibit such 
singular behavior, and how to characterize these singu- 
larities. It is natural to ask whether these singularities 
can appear when no bulk drive is present. 

In this Letter we study boundary-driven diffusive sys- 
tems and show for the first time that singular behavior 
in (j> can indeed occur in these systems. We first provide 
an example of an exactly solvable model which exhibits 
such a singularity. Then we use numerics to demonstrate 
that this phenomenon occurs in a broad class of mod- 
els (for example, in a boundary-driven Ising model) . We 
elucidate the mechanism by which the singularities arise, 
by introducing a simple model which contains the essen- 
tial ingredients leading to the singularity. This model 
shows the close connection between the singularities and 
a spontaneous symmetry breaking. We provide guide- 
lines to systems which can be expected to show this sin- 
gular behavior. Throughout our analysis we use analogies 
with previous works on the Fokker-Planck equation in the 
presence of small noise. Singular behaviors in such mod- 
els have been studied extensively in the past [2ll - |26| , and 
we draw analogies with these works whenever possible. 

Background theory - We study large deviations of bulk- 
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conserving diffusive systems which are driven out of equi- 
librium by the boundaries. For simplicity we consider 
one dimension, with < x < 1. The conserved density 
p (x, t) is related to the current J (x, t) by d t p + d x J = 0, 
where the current is given by 



J 



-D (p {x, t)) 8 x p {x, t) + ^Jo (p {x, t))r, (x, t) . (1) 



D(p(x,t)) is a density-dependent diffusivity func- 
tion, while a{p{x,t)) controls the amplitude of the 
white noise r)(x,t), which satisfies (r)(x,t)} — and 
(r? (x, t) r] (x', t')) = N-i-5 (x - x') 5 (t - t'). The prefac- 
tor A -1 in the noise variance results from the fact that 
we have scaled distances by the system size N, and time 
by A 2 . After this rescaling the noise is small due to the 
system size, as a consequence of the coarse graining and 
rescaling. Eq. ([T]) describes a broad range of transport 
phenomena, including electronic systems, ionic conduc- 
tors, and heat conduction 0,H3,ll3- D (p) and a (p) are 
related via a fluctuation-dissipation relation, which for 
particle systems reads a (p) = 2kBTp 2 K (p) D (p) where 
k(p) is the compressibility [2]. For example, for diffus- 
ing hard-core particles, D = 1 and a = 2p(\ — p) 
The system is attached to reservoirs which fix the den- 
sities at the boundaries of the segment, p(0) = pl and 
p(l) = pr. If pl ^ Pr a current is induced through the 
system, driving it out of equilibrium. 

The average, or most probable density profile for the 
system p, is obtained by solving d x [D (p) d x p] = 0, with 
p(0) = pl and p(l) = pr at the boundaries. In equilib- 
rium (i.e. when pl = Pr), the steady-state probability 
of any other density profile p (x) is easy to obtain - the 
LDF <j) [p] is then given by the free-energy which is local 
in p, (f)[p\ = j f(p, p) dx, where 



f(p,r) 



dpi 



» 2D(p 2 ) 
dpi — -— - 



(2) 



Note that in this case p is constant, p = pl = Pr- By con- 
trast, the steady-state probability distribution away from 
equilibrium is notoriously hard to compute, and very dif- 
ferent from the naive guess <f>\p] = J f (p, p) dx, now with 
space dependent p(x). In fact, as stated above, <f>\p] is 
non-local. 

To compute the large deviation for the model de- 
scribed above, one uses the fact that the probability of 
a history {p(x, t) ,J(x,t)} during time — oo < t < is 
P ~ cxp (-NS), where the action S is given by [J0| 



S = 



f 1 [J(x,t) + D(p(x,t))d x p(x,t)] 2 
Jo X 2a(p(x,t)) 



(3) 

For large N, the probability P ~ e N ^ p fi is given by 
4>[pf] = infp,./ 5* , where the infimum is over histories sat- 
isfying dtp + d x J = 0, with initial and final conditions 
p(x,t — > — oo) = p (x), p(x,t = 0) = pf (x), and bound- 
ary conditions p(x — 0,t) — p L and p (x = = pR, 

In many cases, including in equilibrium and in pre- 
viously studied exactly solvable non-equilibrium models 



0, [Io[, the action S in Eq. (J3j) has a single local mini- 
mum and <j>[pf] is then a smooth functional. However, 
as we show below this need not be the case. The ac- 
tion can, in general, have more than one local minimum 
in the space of histories {p (x,t) , J (x,t)}. In regions 
of the space of final states where pf has more than one 
minimal history leading to it, the global minimum might 
switch between two locally minimizing solutions. This 
is analogous to a first-order phase transition in equilib- 
rium statistical mechanics, where the system switches be- 
tween two metastable states, which are both local min- 
ima of the free-energy. The transition between local min- 
ima is accompanied by a jump in the functional deriva- 
tive of the large-deviation 6<p/8p. We will therefore refer 
to it as a Large Deviation Singularity (LDS). This phe- 
nomenon, first studied by Graham and Tel [21] , is unique 
to non-equilibrium: in equilibrium 4> is smooth whenever 
the dynamical model (Langevin equation) contains only 
smooth functions. Note that while LDSs are expected 
to be generic in models where the zero-noise dynamics 
feature a number of basins, or unstable fixed points, here 
the only fixed point is p(x), and therefore the existence 
of the singularity is not guaranteed, and indeed is not 
present in the previously studied models [1, H, [To| ■ 

We first study this phenomenon in an exactly solvable 
model, and then consider its generalizations. 

Exactly solvable model - Consider the model in Eq. (fTJ) 
with D = 1 and a quadratic a (p) = 1 + p 2 , a parabola 
clear above the axis [29[ . In Q it was shown that the LDF 
is given by 4>[pf] — min0 ea:t , where <p ext are extremal 
values of the action given by 



bext = dx 
Jo 



f(Pf (x),g(x)) - In 



g'(s) 
P (x) 



where / (p, g) is defined in Eq. ([2]) and g (x) is an auxil- 
iary function satisfying the differential equation 



= 



g (x) - p f (x) g" (x) 



°(g (»)) 



W (*)] 



(4) 



with boundary conditions g (0) = pl, and g (1) = pR. 
Note that as D = 1, the most probable configuration 
p(x) is linear, with p(0) = pl and p(l) — pr. As we 
now show, solutions to the differential equation Eq. ([4| 
with initial and final boundary conditions may be non- 
unique, i.e. there exist profiles pj (x) for which more than 
one solution g (x) exists. Eq. (|4]) is solved via a numerical 
shooting method [3(j: here we treat the problem as an 
initial value problem with initial conditions g (0) = pl , 
and g' (0) = c, and scan systematically over values of c to 
find solutions with g (1) = pr. This type of exhaustive 
search ensures that all extremal states are discovered. 

To illustrate the existence of multiple solutions, 
we consider profiles of the form Pf (x) — p{x) + 
cui cos (ttx/2) + a2 sin (ttx), varying ot\ and and with 
boundary-conditions pl = —3, pr = 3. Fig. []Ja) shows 
an example of a profile pf for which only one solution 
g (x) to Eq. Q exists. In contrast, in Fig. 0Jb) a profile 
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FIG. 1. The model D = 1, a = p 2 + 1. (a) A profile p f (x) 
(solid line) with one extremal solution <f> [p]. The correspond- 
ing g (x) also plotted (dashed line), (b) A profile pf (x) with 
three extremal solutions, (c) Number of extremal solutions as 
a function of pj (1/3) , pf (2/3). Gray region: three solutions. 
White region: one solution. Also shown are the caustics (solid 
line), and the switching line (dashed line). 



p f (x) with three solutions g (x) is shown, two of them 
corresponding to local minima, and one to a saddle point. 
Fig. [BJc) shows the region in which there are multiple so- 
lutions. Here we have chosen to parametrize the profiles 
Pf (x) in terms of pf (1/3) and pf (2/3), which are simply 
related to a± and ct2- Note the marked caustics, indicat- 
ing the boundaries between regions with one and three 
extremal solutions, and the switching line, on which the 
two locally minimal solutions have the same value for the 
action. On the switching line the gradient of the LDF is 
discontinuous 2l|; the occurrence of this LDS is the fo- 
In addition, the history preceding a 



cus of the paper, 
rare event is expected to be different on both sides of 
the switching line, as the history minimizing the action 
changes. The cusp is found for profiles pf relatively "far" 
from p, and for pf (1/3) positive and pj (2/3) negative 
(here pl < Pr)- More generally, the phase space of pro- 
files pf is infinite dimensional, the caustics and switching 
line become manifolds. The picture shown in Fig. [TJc) 
is a particular two dimensional cross section. We find 
similar behavior for p/ of similar shapes which are not of 
the exact form described above. 

In fact, for this model one can prove that: (a) for any 
non-equilibrium boundary conditions (pl ^ Pr) there 
exist profiles pf for which there is more than one solu- 
tion to Eq. (01, and (b) profiles which have two locally 
minimizing histories with the same value of <p ex t always 
exist. A proof of these facts will be given elsewhere [3l| . 

The existence of multiple minimizers of S can be intu- 
itively understood as follows: looking at Eq. (|3]), we see 
that the contribution to the action is smaller wherever 
a (p) is high. If the variation in a (p) is large enough, 
trajectories passing through regions of high a (p) may be 
favored. If there are two such regions, as around a mini- 
mum of (j, there might be different paths of the action uti- 
lizing the different favorable a (p) regimes. When D (p) 
also varies, we expect the same logic to apply to regions 



with high and low a (p) / D (p) . This argument suggests 
that the phenomenon is robust, and will occur in other 
modes with similar features, i.e. when a (p) /D(p) has 
a pronounced minimum. Below we make this argument 
precise, but first we demonstrate the generality of the 
phenomenon by studying it on a different model which 
admits a concrete microscopic realization. 

Boundary driven Ising model - We turn to study a 
boundary-driven Ising model, a lattice gas with on-site 
exclusion and nearest-neighbor interaction. (This is a 
variant of the Katz-Lebowitz-Spohn model [32j , but with 
no bulk bias.) Each site i = 1, .., N of a one-dimensional 
lattice can be either occupied ("1") or empty ("0"). The 
jump rate from site i to site i + 1 depends on the occupa- 
tion at sites i — 1 to i + 2 according to the following rules: 

0100 ^ 0010,1101 ^ 1011,1100 y> £ 1010,1010 ^ 
1100, and their spatially inverted counterparts with iden- 
tical rates. The parameter < e < 1 corresponds to at- 
tractive interactions between the particles; S controls the 
density dependence of the mobility. As shown in [H, HH , 
for each parameter set (e, 5) one can write implicit an- 
alytic equations for D (p) which can then be inverted 
numerically. Then a (p) is obtained via the fluctuation- 
dissipation relation. Fig. [Ua) shows D (p) and a (p) for 
(e, 5) = (0.05, 0.995). For equilibrium boundary condi- 
tions this model admits an Ising measure. 

We now solve for local minimizers of the action S, us- 
ing the numerical technique described in [11| . The nu- 



merical solutions are obtained by gradually changing the 
end profile pf, while continuously maintaining a locally- 
minimizing history p(x,t). Different locally-minimizing 
solutions are obtained by changing the final profile to 
enter the bi-stable region from different directions, see 
Fig. [2]^c). As for the exactly-solvable model, we again 
find configurations p/ (x) for which multiple local min- 
imizers of S exist. In Fig. HJb), we show two different 
histories p(x,t) leading to the same profile p/, which is 
chosen again to be of the form: pf = p (x) + a± sin7rx + 
cti sin 2-7TX. In Fig. (2fc) we plot the numerically-obtained 
region in the p/(l/3),p/ (2/3) plane with multiple min- 
imizers, together with the "switching line" and the caus- 
tics. In addition, for a particular pf (x) we depict the two 
histories leading to it by showing the time-dependent val- 
ues of p (x = 2/3, t) against p (x = 1/3, i). Lines of equal 
LDF are plotted in Fig. [2Jd) to show the jump in its 
gradient. 

Note that in this model and for the chosen parameters, 
a(p) /D (p) has a double- hump structure with a deep 
minimum, so an LDS is expected from the simple consid- 
erations discussed above. Indeed, we have experimented 
with various forms of D (p) and a (p), and conclude that 
an LDS occurs when a (p) / D (p) has a minimum which is 
deep enough. Recall that by fluctuation-dissipation is re- 
lated to the compressibility a (p) /D (p) = 2ksTp 2 K (p). 
These features have to be large enough for this to hap- 
pen; small changes to a model which does not feature an 
LDS are generally not enough. When LDSs do appear, 
we have always found them in regions of phase space 
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FIG. 2. The boundary-driven Ising model, (a) D (p) and 
<j(p). (b) Evolution of two locally-minimizing histories lead- 
ing to the same pf. (c) The coexistence region, showing the 
caustics (solid black line), switching line (dashed line), and 
cross-sections of two histories (gray line). Dashed arrows de- 
pict paths of the final state in the numerics which yield differ- 
ent local minimizers. (d) Contours of equal LDF (solid lines) 
and the switching line (dashed). 



FIG. 3. (a) Passing a mass through the noise barrier p £ 
[—e/2, e/2]. (b) Two profiles which can be reached by only 
pushing mass "downhill" through the noise barrier, (c) A 
symmetric history leading to pf requires passing a mass 2Af , 
shown in gray, "uphill", (d) A symmetry-breaking history 
which requires passing only a mass M "uphill". In (b-d) the 
straight thin line depicts p, and the bold lines pf. 



where the profiles p f (x) have a similar shape to that 
shown in Fig. HJb) and Fig. [2][b). The exact conditions 
and profiles for an LDS to appear must be studied on a 
model-specific basis. Below we provide a simple model 
which illustrates the mechanism leading to LDSs in mod- 
els with these features. 

Mechanism - Finally we elucidate the connection be- 
tween a deep minimum in a (p) and the LDS in the con- 
figurations discussed above. To do so we introduce a 
simple sweater sleeve model. We consider a model with 
D = 1 and a (p) of order one everywhere except for a 
narrow range of density values, where it is small: 



a(p) = 



P e [-e/2, e/2] 

otherwise 



Here e is a small parameter, and o\ (p) is some func- 
tion independent of s. Then a (p) has a deep minimum 
around p = 0. (The fact that a (p) is discontinuous is 
not essential - a smoothed version can also be used.) 

The key to estimating the action is noting that pushing 
a mass through the density region p e [—e/2, e/2] may 
be costly in terms of the action, creating a "noise bar- 
rier" . Consider the cost of passing a small mass element 
m from p = —e/2 to p — e/2 or vice versa, by applying 
a current J, see Fig. Ela). This process is done over a 
time At. The region of space where p £ [—e/2, e/2] is 
of length Ax. As J = m/At,d x p = e/Ax, the action 
S = J dxdt (J + d x pf / (2a) is, to order O (e^ 1 ), 



AxAt 



(- 



2a e 2 \At Ax 



2<r £ 2 



m v H h 2me 



where v = . We distinguish between two cases: J 
"uphill" (against Fick's law, requiring a strong noise), i.e. 



sign(J) = sign (d x p), and "downhill" with sign(J) = 
—sign(d x p). In the "uphill" case, and have the 
same sign, and minimizing 5" over v we obtain v = e/m 
or S > 2si -f O (e°) . To push a macroscopic mass M , 
the bound will then read S > — + O (e°) . In contrast, 
if J is "downhill" , J + d x p can be made small, with no 
bound of order O (e -1 )- 

We now argue that for some profiles pf, such as the one 
depicted in Fig. [3jc,d), the global minimizing history is 
not unique. We consider boundary conditions pl < —e/2 
and e/2 < pa, so that pf crosses the noise barrier three 
times. To highlight the symmetry-breaking aspect of the 
phenomenon, we focus on a model with a Z2 symmetry 
where o\ (p) = o\ (— p), p^ — —pn, and on a final pro- 
file satisfying pf (x) = —pf (1 — x). Then a (p) ,D and 
Pf have a symmetry under the combined operation of 
p — > —p and x — > 1 — x. Let p svm {x, t) be the minimizer 
of S subject to this symmetry. In the space of histories it 
is extremal, but not necessarily the minimizing history. 
Indeed, we show that other solutions have lower action, 
spontaneously breaking the symmetry. Referring to Fig. 
[3fc) we note that in a symmetric history p(x,t) we have 
p(x = 1/2, t) = 0, and masses in the region shaded in 
gray, 2M, must by pushed "uphill" through the p = 
line, hence S > However, in the symmetry-breaking 
history shown in Fig. EJd) only a mass M has to be 
pushed "uphill", so that S = fjJf + O (e°). Therefore, for 
a small enough e this solution is favorable. By symmetry, 
the solution with p — > —p and x — > 1 — x has the same 
action. Therefore the symmetry is spontaneously broken, 
leading to the LDS with a switching line on profiles obey- 
ing this symmetry. Note that some profiles pf, such as 
those in Fig. EJb), can be generated by a history which 
only pushes mass "downhill" through p £ [—e/2, e/2], 
and the above argument for multiple histories will not 
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hold. While the above argument emphasizes the break- 
ing of a symmetry, the phenomenon exists even in the 
absence of an explicit symmetry in either the final con- 
figuration or the model, as was demonstrated in the first 
two models considered. 

The non-differentiability shown in Figs. [IJc) and 
HKc,d) is perhaps the simplest possible singularity [26| . 



Due to the high dimensionality of the phase space of 
profiles, more elaborate structures are possible. Under- 
standing and classifying the possible singularities in these 
models is an exciting direction for future research. 
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